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This paper presents the work devoted to the study of the operation of a miniaturized 
membrane Stirling engine. Indeed, such an engine relies on the dynamic coupling of the 
motion of two membranes to achieve a prime mover Stirling thermodynamic cycle. 
The modelling of the system introduces the large vibration amplitudes of the membrane 
as well as the nonlinear dissipative effects associated to the fluid flow within the engine. 
The nonlinearities are expressed as polynomial functions with quadratic and cubic 
terms. This paper displays the stability analysis to predict the starting of the engine and 
the instability problem which leads to the steady-state behaviour. The centre 
manifold-normal form theory is used to obtain the simplest expression for the limit 
cycle amplitudes. The approach allows the reduction of the number of equations of the 
original system in order to obtain a simplified system, without loosing the dynamics of 
the original system as well as the contributions of nonlinear terms. The model intends to 
be used as a semi-analytical design tool for the optimization of miniaturized Stirling 
machines from the starting to the steady operation. 

© 2009 Elsevier Ltd. All rights reserved. 


1. Introduction 

One of the new application of the Stirling engine is the free piston configuration proposed by Beale of the Ohio 
University by the end of 1960s [1 ]. Without any mechanical linkages between the piston and the displacer, this architecture 
allows the realization of fully sealed machine which leads to the most reliable machines compared to the more classical 
mechanical arrangements. In such a machine, the reciprocating motions of the piston and displacer are effected by the 
dynamical and the fluidic interactions. The operation of free piston Stirling engine (FPSE) results from a self starting 
characteristic followed by a stable operation stage when nonlinear phenomena modify the initial linear behaviour limited 
to the small amplitudes of the moving parts. These nonlinearities are generally related to pressure drops when the working 
fluid goes back and forth within the heat exchangers and the regenerator and consequently allow the stable operation. 
Thus, the thermodynamic cycle can be achieved with appropriate parameters. The main advantages are: a simple 
mechanical design; minimal lateral forces on moving components which minimizes the engine wear; high energetic 
efficiency; great life time and low cost manufacture. 

Without any mechanical linkage the thermodynamic optimization is a complex issue. Principal design parameters, 
swept and dead volume ratios, phase angle and speed of engine are coupled. In depth dynamical analysis is therefore 
necessary to account for design parameters effects on performances. Though a numerical model is to be considered usually, 
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a semi-analytical approach would allow an efficient preliminary design for a strongly coupled problem with numerous 
design parameters [2]. The design of such an engine requires to master complex thermal, mechanical and fluid dynamics 
interactions. In most analytical studies assumptions of the isothermal Schmidt cycle [3] are adopted then the traditional 
approach consists in linearizing of each force which acts on the moving parts and considering the whole problem as a 
coupled two degrees of freedom (DOF) system. One obtains therefore the dynamical equations in the neighbourhood 
of a stable operating point which leads to the dynamical and consequently the thermodynamic characteristics of 
the FPSE [4-6]. 

Few works deal with miniaturized Stirling machine for which the moving components are replaced by vibrating 
membranes. Bowman in [7,8] developed an analytical model to design membranes which match the dynamical criteria 
(notably swept volume and intrinsic dynamical characteristics). However, this analysis is restricted to linear behaviour 
whereas the geometric nonlinear phenomena appear very soon as a flexural behaviour is to be considered [9]. Under these 
conditions, the approach proposed by De Monte Refs. [4,5] and Rogdakis [10] which includes the nonlinear dissipation 
phenomena associated to the fluid flow is completed hereafter by the mechanical nonlinearity of each membrane. The later 
appears to be a crucial point to forecast the performances of a membrane Stirling engine. 

During recent years, lots of works which aim at the understanding of the dynamic behaviour of systems with nonlinear 
phenomena have been developed. Perturbation methods, such as averaging or the methods of multiple scales have been 
used in many studies for simplification purpose of the effects of nonlinear phenomena which occur in many problems 
related to engineering and science [11 ]. The centre manifold approach is the one which leads to the most simplified model. 
It accounts for the local bifurcation behaviour in the neighbourhood of a fixed point of the nonlinear system [12]. 
The centre manifold approach reduces the number of equations of the original system in order to obtain a simplified 
system without loosing the dynamics of the original system as well as the effects of nonlinear terms [13]. Moreover, the 
normal form approach can be jointly used to eliminate as many nonlinear terms as possible through nonlinear successive 
changes of variable. 

The present paper depicts the application of the centre manifold theory and the normal form approach to a membrane 
miniaturized Stirling engine (MMSE) under the Schmidt assumptions. The membrane micro-engine arrangement and its 
expected behaviour are first introduced. Next, the dynamic of the membranes will be outlined. A first simplification 
method based on a modal projection leads to a set of two coupled second-order differential equations which is 
representative of the system’s dynamic. The coupling from the isothermal evolution assumption which accounts for the 
pressure to be analytically related to temperature and volume variations is established. The system stability is then 
assessed by the linearization of the equations at a fixed point and self starting conditions are outlined. Finally, the centre 
manifold and the normal form approaches will be used in order to predict limit cycle amplitudes and the estimation the 
performances of the engine eventually. 

2. Mechanical arrangement and qualitative working of a membrane miniature Stirling engine 

Based on the patent of a Stirling micro-machine from Bowman [7], Fig. 1 shows its generic mechanical arrangement as 
well as the relevant parameters for the study. 

The operation of the engine is based on the bending of two membranes called m exp and m comp which can therefore 
modify the volumes of the chambers V ei V t and V c , \/ s , respectively. The compression membrane denoted m comp ensures the 
role of the piston and the expansion membrane denoted m exp the role of the displacer of the Stirling machine. For an 
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Fig. 1 . Mechanical arrangement of the micro-Stirling engine. 


























































































796 


F. Formosa / Journal of Sound and Vibration 326 (2009) 794-808 


Table 1 

Geometric characteristics used for this study. 
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Fig. 2. Representative stages of operation. 


engine, the volume V e is at the highest temperature Th whereas V c and V* are maintained to lower temperatures Tc and Ti 
thanks to the cold exchanger. On can note that we make a distinction between the spaces above and under the cold heat 
exchanger (with respect to the z-axis shown in Fig. 1 ) in order to account for a pressure drop as well as a temperature 
difference between the two. A moving regenerator V Reg provides thermal insulation between V e and V,- and allows better 
performances alternatively storing and releasing heat to the working fluid. Table 1 gives the geometric characteristics used 
for this study. 

The representative stages for the operation of the engine can be described from a reference position of the bending 
components. In Fig. 2a, the compression membrane is at its lowest position as the expansion membrane is at its highest so 
both V c and V t are maximum and V e is minimal. In Fig. 2b, the upward motion of the compression membrane induces a 
global compression of the working fluid. The subsequent temperature increase is balanced by an outward heat flow at the 
cold exchanger. With a downward motion of the expansion membrane the working fluid is partly forced trough 
the regenerator. Therefore, heat energy is provided to the fluid by the regenerator material so it ideally emerges from it at 
the temperature Th (Fig. 2c). The global pressure consequently rises as for a constant volume the temperature becomes 
higher. Submitted to this pressure, the compression membrane moves to its initial position and by doing this, some 
mechanical power is extracted from the engine. Spring forces lead to an out of phase motion of the expansion membrane 
(Fig. 2d). The heat source provides the amount of energy to ideally keep a constant high temperature Th during 
the resulting global expansion. At the end of the expansion stage the expansion membrane is at its initial position and the 
regenerator material recovers the amount of heat energy harnessed by the fluid during the previous blow (from Vi to Ve). 
In the ideal case, the fluid went back to the compression spaces (V c and V,) at the temperature Tc. 


3. Semi-analytical model of the membranes 

In order to simulate the global behaviour of the engine, one first interests in the modelling of the moving parts. For the 
considered micro-engine the compression membrane is assumed to be a single thin, circular and clamped plate whereas 
the expansion membrane is supposed to be composed of two clamped circular shells linked by a rigid stem at their centre 
(see Figs. 1 and 13). 

Haterboucha-Benamar [9,14] developed a theoretical model based on Hamilton’s principle and spectral analyse. A semi- 
analytical relation between the eigenmodes, the eigenfrequencies and the amplitude of bending for large deflections can 
thus be established. This method has been successfully applied for numerous typologies of structures such as simply 
supported or clamped beams and plates [15-17]. The method is especially efficient not only to account for the 
amplitude-frequency relation but also represents the effect of large deflexion on the shape of the membrane. Because of 
the importance of the swept volume on the performance of the micro-engine, this last result happens to be relevant in the 
design procedure. The main assumptions used for the modelling of the membrane are the following: 


• Membranes are submitted to mechanical solicitations such as their strains are only axisymmetric. 

• The only displacement under consideration is the displacement in the z direction denoted w*(r). 

• Plane stresses are considered (shear stresses are neglected). 

• Rotational kinetic energy is neglected. 
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Fig. 3. Normalized resonance frequency with respect to normalized amplitude. 


To the extent that the harmonic motion is supposed, strain and kinetic energy can be obtained using the basis of linear 
eigenmodes which form the spectral basis to express the transverse displacement. In the case when large deflexion is to be 
considered, nonlinear geometric behaviour leads to a tensor of the fourth order associated to in plane induced tension in 
addition to the classical stiffness and mass ones. It follows from the application of the Hamilton’s principle the sets of 
algebraic nonlinear equations which can be numerically solved. Provided that the spectral basis can be written using 
Bessel’s functions of the first and second kind, this approach has the advantage to provide a semi-analytical model of the 
nonlinear dynamical behaviour for a membrane. 

In the following, the method is applied to the compression and the expansion membranes although in this study, 
the restriction of the spectral basis to a single vector will be used in order to obtain only one DOF equation for each of the 
membrane such as 


<h m !i + <7i fc ii =p/? (1) 

In which cjq and qi are the acceleration and displacement response of the generalized coordinate and p is a pressure applied 
to the membrane area. m\ v k\^ y andare the mass, linear stiffness, nonlinear stiffness and reduced pressure 

coefficient terms which are evaluated from the spectral analysis of the membrane equilibrium equation (the method is 
described in Appendix A). 

It is worthy of note that the previous Eq. (1) can be seen as a conservative Duffing’s oscillator from which the 
characteristic evolution of the resonance frequency of the membrane with respect to its deflexion is plotted in Fig. 3. 

Yet, the single DOF model is relevant to represent the effect of the amplitude of deflexion on the first resonance 
frequency for a membrane under large strain assumption. Thorough study with larger spectral basis would not involve 
much theoretical difficulties but would consequently increase the number of DOF. 

The development of the analytical method for the compression membrane and the expansion membrane is developed in 
Appendix A. 

4. Analytical model of the MMSE 


The governing equations of the micro-engine can thus be obtained using a two steps procedure. The first step as was 
shown before establishes a single DOF nonlinear equation for each of the membrane such that 
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with the following initial conditions: 



q = 0, b = 0 and q, = 0, q . =0 (4) 

~P\t=0 —P\t=0 -d|t=0 -d|t=0 

where (m p , k p , b p ) and (m d , k df b d ) are the mass, linear stiffness and nonlinear stiffness terms which are evaluated from the 
spectral analysis of the compression membrane and the expansion membrane, respectively, pe, pi, pc and ps are the gas 
pressure of the chambers as defined in Fig. 1. 

Eq. (2) stands for the dynamic behaviour of the compression membrane. q p is the displacement response of the 
compression membrane DOF and q p is the derivative of q p with respect to time. £ p is a modal damping coefficient added to 
represent the mechanical dissipative phenomena. It is worthy of note that the damping coefficient £ p can represent the 
effects of a linear load applied to the micro-engine. As an example, such a load can be induced by an electromagnetic 
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transducer. As a consequence, its value will be chosen in a quite wide range. Set of Eqs. (2)-(4) can also be seen as two 
oscillators coupled by mechanical forces. 

Eq. (3) stands for the dynamic behaviour of the expansion membrane. q d , q d and q d are the acceleration, the velocity 
and displacement response of the expansion membrane DOF, respectively. £ d is a modal damping coefficient added to 
represent the mechanical dissipative phenomena. 

For the second step, under usual assumptions of isothermal evolutions and ideal working fluid, the instantaneous 
pressures within the chambers can be analytically expressed. The total mass of gas inside the various chambers is defined 
as the sum of the mass of the gas contained in each of the chamber defined in Fig. 1 : 


M = m e + m h + m Reg + m* + m k + m c (5) 

If we substitute the expressions of the masses in Eq. (5) with the ideal gas law (e.g. m e = pV e /RTh), we obtain 


p = MR 



Ve V h VReg Vi V k V c 

Th Th Tr Ti Tc Tc 



where R is the considered gas constant per unit of mass and Tr the mean temperature of the regenerator. As expected, the 
pressure depends on the motion of both membranes. As an example, the motion of the compression membrane modifies 
the pressure. Moreover, the different pressures of the chambers are linked by the pressure losses effects which occur in the 
heat exchangers as well as in the regenerator. 

Thus denoting p the mean global pressure, we get: 


pe = P + ^ Ap Reg , pc = p -1 Ap Ech and pi = p - i Ap Reg +i Ap Ech . 
where pe, pi and pc are the pressures defined in Fig. 1. 

The expressions of the pressure losses Ap Reg and Ap Ech are highly nonlinear relation in which the classical Reynolds 
number is the main variable. To simplify the approach, cubical nonlinearities are deduced from the classical pressure loss 
expressions [18,19] by Taylor’s development. Therefore, the general form of Eqs. (2) and (3) can be given in the following way: 


~Qd~ 


■ 0 

1 

0 

0 ■ 


~Qd~ 


0 

Qd 


~ S dd 

~ D dd 

~$d P 

~^d P 


Qd 

+ 

fNL d (q d , q p ) 

Q P 

— 

0 

0 

0 

1 


Q P 

0 

Qp 


~S p d 

~D p d 

-S P p 

-D pp _ 


Q P 


fNL p (q d , q p ) 



with the initial conditions: q P]t _ Q = 0, q P]t _ Q = 0 and = 0, q^ t=Q = 0. Wherein, Sy and Dy are the stiffness and the 
damping coefficients between the i and the j DOF. fNL d (x,y) and fNL p (x,y ) contains the nonlinear stiffness terms. 
Hereafter, Eq. (7) will be written in state variable x l = [q d ; q d ; q p ; q p ]: 


where 
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, b d]), x<8>x defines the basis of quadratic terms 2 and x<8>x<8>x 


5. Solution procedure 


The arising of steady oscillations in coupled nonlinear oscillators is classically associated to a supercritical Hopf 
bifurcation. At the fixed point also called the bifurcation point and denoted x 0 hereafter a couple of complex conjugated 
eigenvalues associated to the system crosses the imaginary axes. Hence, this point becomes unstable and because of the 
nonlinearities of the system a stable limit cycle can occurs in its vicinity. 

5.2. Unstable fixed point 


The fixed point x 0 is obtained by solving the nonlinear static equations for a given set of parameter. It follows from 
Eq. (8) provided x = 0 that the only fixed point is x^ = [q d0 ; 0; q p0 ; 0] = [0; 0; 0]. 

The stability is investigated by calculating the eigenvalues of the constant Jacobian matrix at the fixed point x 0 . 
It is possible to obtain the fourth-degree characteristic polynomial for which the roots are these eigenvalues, we set 

r 4 + ar 3 + ffr 2 + yr + 3 = 0 


(9) 
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where 


a = D pp + D dd 

P = DppDdd ~ D pd D dp + S dd + Spp 

y = ~^dp$pd ~ Dpd^dp "T ^dd^PP "T ^PP^dd 

<5 = SppS dd — Sp d S dp (10) 

Hence, depending on the operating parameters, three situations can occur: 

If at least one of the real parts of the four roots would be greater than zero the system would be unstable in the vicinity 
of the fixed point. 

If every real part would be less than zero the system is asymptotically stable. The motions of the membrane will tend to 
zero within a finite period of time. 

If two roots have zero real part and the other two have negative real part, steady oscillations can occur. 

The operation of the MMSE relies on the last case for which the roots of Eq. (9) at the fixed point can be expressed as two 
complex conjugates expressions: 

%ci = -jVv7a = -j®i; 4c2 = j VyJv = +j®i (ii) 

Ajsi = l/2(-a - ji^-a 2 +4<5a/y) = a 2 -ja» 2 ; X }s2 = l/2(-a + j^/-a 2 +4<5a/y) = a 2 +jffl 2 (12) 

where j is the complex number such as j = ^/^T. 

Moreover, applying the Routh-Hurwitz criterion applied to the characteristic Eq. (9) results in the three following 
conditions for the stability: 

Cl : a>0 (13) 


C2 : a/? — y >0 (14) 



These conditions are well known for the analyses of the FPSE [6]. De Monte et al. established the so called P4 theorem [4,5]. 
The evolutions of the terms of Eqs. (13) and (14) with respect to Th/Tc ratio and p/p 0 ratio where p 0 is a reference pressure 
are plotted on Fig. 4. 

As can be seen on the two graphics of the left, the first two conditions are usually met. The third condition defines an 
unstable region related to the highest temperature and mean pressure parameters. In the case of two pairs of complex 
conjugate eigenvalues, the evolution of their real parts is studied for increasing values of Th. The determining rule of the 
temperature for the apparition of the loss of stability is underlined in Fig. 5. As a result, a critical starting temperature can 
be defined. 

The effect of the load of the system represented by the damping coefficient is reflected in Fig. 6 in which increasing £ p 

leads to the stabilization of the system. Thus, the load may have to be controlled during the starting phase to make it easier. 

5.2. Centre manifold approach 

The linear analysis cannot represent the behaviour of the engine after the starting point. Therefore, another theory has 
to be used to study the post bifurcation behaviour. This section describes the application of the centre manifold [13] 



Fig. 4. Stability criteria for the behaviour of the engine. Left: Cl: a>0, middle: C2: a/?-y>0 and right: C3: (y/a) 2 —/3(y/oc)+<5>0. 
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(a) Im (b) Im 




Fig. 5. Evolution of the system’s poles with respect to temperature ratio Th/Tc: (a) first positive imaginary part pole and (b) first positive imaginary part 
pole. 


(a) . (b) 




Fig. 6. Evolution of the system’s poles with respect to £, p : (a) first positive imaginary part pole and (b) first positive imaginary part pole. 


method to simplify the system defined by Eq. (8). Locally, the stability of the centre manifold is equivalent to the stability of 
the original system. Again, the bifurcation appears when one or several eigenvalues cross the imaginary axis in the complex 
plane with the variation of a control parameter (Figs. 5 and 6). 

At the Hopf bifurcation point, one defines the coordinate transformation matrix T such as 
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where ai, co i and a 2 , co 2 are the real and imaginary parts of the complex eigenvalues defined in Eqs. (11) and (12). Note that 
at the bifurcation point ai = 0. Thus, the previous system Eq. (8) can be written again in the form: 


Vc = J c Vc + G 2 (V C , Vs) + G 3 (v c , Vs) 


V s = J S V S + H 2 (V C , V S ) + H 3 (V C , V s ) (17) 

where J c and J s have eigenvalues such as Re(2j 
components are polynomials of degrees 2 and 3 in the components of vectors v c and v s namely v cl , v c2 and v s i, v s2 . 

By considering the physically interesting case, it may be assumed that a 2 is negative. v s defines a stable subspace which 
is the linear approximation to the stable manifold. When we choose an initial condition on this stable manifold sufficiently 
close to the fixed point the solution curve will go toward the fixed point. In accordance to the parameters, an unstable 
subspace defined by the eigenvalues of A such as their real part is zero can exist. The centre manifold theorem allows to 
represent locally the centre manifold as {[v c , v s ] such that v s = h(v c ), h(0) = 0, Dh(0) = 0} and consequently reduce the 
initial problem to a two dimensional one. 

Substituting into the second line of Eq. (17) one obtains 

v c =J c v c + G 2 (v c ,v s ) + G 3 (v c ,v s ) ( 18 ) 


) = 0 and Re(2j s ) = a 2 <0. G 2 , G 3 , H 2 and H 3 are vector for which the two 


D vc (h(v c ))v c = J s h(v c ) + H 2 (v c , h(v c )) + H 3 (v c , h(v c )) 


(19) 
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Eq. (19) cannot be solved analytically. However, it is possible to define an approximate solution of h by a power expansion 
without constant and linear terms. Thus, h is expressed as a power series in the components of v c of degree m: 

m p p . 

h< v c) = E ( 2 °) 

P—i+j—2 j =0 1=0 

Replacing v c in Eq. (19) one obtains Eq. (21) form which the resolution provides the complex coefficients a,j of the terms of h. 

D vc (h(v c ))(J c v c + G 2 (v c , h(v c )) + G 3 (v c , h(v c ))) = J s h(v c ) + H 2 (v c , h(v c )) + H 3 (v c , h(v c )) (21) 

As a result, the dynamic behaviour of the system is determined by a reduced two DOF system given by 

Vc = J c Vc + G 2 (V C , h(v c )) + G 3 (v c , h(v c )) (22) 


5.3. Normal form analysis 

The normal form analysis aims at a transformation of a system of nonlinear equations through a sequence of nonlinear- 
near identity transformations to eliminate as many nonlinear terms as possible. Those which cannot be removed are called 
the secular or the resonant terms. This simplest form of the equations is called the “normal form”. The analysis of the 
dynamics of the normal forms yields a qualitative picture of the flows of each bifurcation type. 

The centre manifold theorem has been applied to the initial system and henceforth we consider the resulting Eq. (22) 
detailed below with T c the complex conjugate of v c : 

1>C = Ot-[Vc — CO} Vc + G 2 _\ (Vc, Vc, h(Vc, f'c+))G 3 _l (Vc, Vc^(Vc, Vc)) 


Vc — V C + Ot\Vc + G 2 _ 2 (V C , Vc, h(Vc, Vc)) + G 3 _ 2 (Vc, Vch(Vc, Vc)) (23) 

In which G 2 _i is the first component of the vector function G 2 , where, for k = 0 to 3 and p = 0 to 3: 

G 2 -\(Vc, Vc, h(V c , Vc)) + G 3 _ A (V c , V c , h(V c , V c )) = E v'y ( 24 ) 

2 < k+p < 3 


G2-2(Vc,VcMVc,Vc)) + C 3 _ 2 (Vc,VcMVc,Vc))= E Ppk^c 

2 < k+p < 3 


For the sake of clarity, Eq. (23) can be written as 


(25) 


Vc = Ot\V c — CD 1 Vc +f(Vc, Vc) 


Vc = 00 A Vc + Vc + g(Vc, Vc) (26) 

We observe that if we set z = vc - jv c andz = v c + jv c the previous Eq. (26) is a pair of complex conjugate equations. 
Therefore, one single equation needs to be studied. A Hopf bifurcation is recognized here and the result of these successive 
transforms is well known [13]: 


and 


v c = a^v c - a>i v c + (av c - kv c )(v 2 c + v 2 c ) 


(27) 


16a — (fxxx + fxyy + Sxxy + Syyy) + ^^ I fxy(fxx + fyy) §xy(§xx + §yy) fxxSxx + fyygyy] 

16 b = (fxxy + fyyy ~ Sxyy ~ §xxx) + [2 (f xx + fxy) + ^(fyy ^Sxx) + 5/ xx(f yy ~ Sxy) + 2 g X y ~ fyySxy 

+ 5§XX§yy + 2gyy ~fxy(§XX + 5gyy)] 


(28) 


(29) 


In which k x is the partial derivative to the first order of a function k with respect to the variable x and k xx is the partial 
derivative to the second order of a function k with respect to the variable x. 

Finally, we have 


v c = a- l v c - coi v c + (av c - bv c )(v 2 + v 2 ) 

Vc = CO\ V c + GC-i Vc + (bv c + av c )(v 2 + v 2 ) 

We can transform Eq. (30) in polar form (v c = r cos 6 and v c = r sin 6). Doing so, we obtain 

r = ot^r + ar 3 


(30) 


6 = 0Q\ + hr" 


(31) 
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Fig. 7. Comparison of the results from straight resolution, centre manifold and normal form: (a) expansion membrane generalized coordinate position 
and velocity and (b) compression membrane generalized coordinate position and velocity. 


The amplitude of the limit cycle can be assed by the analysis of the first of the previous equation. The non trivial solution 
establishes the existence condition for the limit cycle as a 1 /a<0. 

The stability of each fixed point can be studied by the sign of the Jacobian J(r). From Eq. (31) J(r) = ai+3ar 2 , therefore 
y(ryi) = ai and J(r/ 2 ) = —2ot\. Thus, if ai is lower than zero and a greater than zero, the only stable fixed point is r f i = 0 and if 
ai >0 and a< 0, the unique stable fixed point is and a limit cycle occurs. 

These two last conditions added to the three previous ones ensure the existence of a stable limit cycle: 

C4:aj>0 and C5 : a<0 (32) 

Phase portraits of Fig. 7 shows the motion of the generalized coordinate for the expansion membrane (q^,^), in Fig. 7a and 
for the compression membrane ( q p ,q p ) in Fig. 7b. The results of a numerical resolution of Eq. 7, the resolution of Eq. (22) 
obtained from the centre manifold approach and the simplest expression of the result from the normal form analysis in 
Eq. (31) are shown in dotted (st), dashed-dotted (vc) and plain line (ap), respectively, and present a good correlation for 
parameters close to the bifurcation point. 

It is worth to notice that only the centre manifold approach followed by the normal form analysis allows an analytical 
expression of the solution. 


6. Results and discussions 

In the previous section the results from the simplification method have been compared to the reference ones from a 
numerical resolution of Eq. (7). As a consequence, the approach is validated and is an efficient tool to study the effects of 
parameters on the motions of the membranes and eventually on the performances of the engine. 

The classical control or variable parameters of a FPSE are: the heat source temperature; the heat sink temperature; the 
load and the internal pressure. The latter is barely controlled as it would require complex devices. In the specific case of a 
membrane Stirling engine when electric generation is considered, the mechanical load on the compression membrane is 
related to an electromechanical converter device. Hence, we assume that it could be a part of the control parameters. 
The heat source temperature can be a control parameter for a quite large range of values whereas the heat sink temperature 
is related to external conditions, temperature specifically. Consequently in the following, we will focus on the effects of load 
and heat source variations. 

For each simulation, the Hopf bifurcation point and the value of the critical ratio of the hot to cold temperature Th/Tc are 
computed. The study is provided near the bifurcation point by choosing the parameter value Th not greater than 1 percent 
of its critical value. Table 2 gives thes values used for the simulations. It is obvious that a lower value than the critical one 
does not match the condition of Eqs. (13)—(15). 

It is observed that the level amplitude is a complex problem. The evolution of limit cycle amplitude is not linear with 
the evolution of a specific parameter. It turns out that the increasing of the load leads to a reduction of the amplitude of the 
compression membrane whereas the expansion membrane amplitude increases. This is further reflected in Figs. 8-10. 

The evaluations of the efficiency and the output power of the system are given in Table 3. As the upper temperature Th is 
varying, the relative efficiency defined as the ratio of the efficiency with the Carnot efficiency q c = l-Tc/TTi, allows a fair 
comparison between the different cases studied. 

It is worth of noting that the ability of the engine to operate as a prime mover relies on the phase between the 
compression and the expansion membrane motions. If the damping coefficient £ p is smaller than 0.08, the compression 
membrane does not lag behind the expansion any more thus leading to a reverse behaviour of the Stirling cycle which is a 
cooler or heat pump. The sign of the mean output power changes in this case. 

From Fig. 11, it can be seen that the efficiency of the engine is strongly dependant on the load. From the small amplitude 
of the compression membrane a smaller compression ratio can be inferred but higher operating frequencies involve higher 
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Table 2 

Evolution of the limit cycles near the bifurcation point with respect to £ p . 



0.04 

0.06 

0.08 

0.1 

0.12 

0.14 

TcritifTc 

1.09748 

1.67958 

2.08963 

2.41272 

2.68125 

2.91167 

Thtest = 1.005 *Thcrit 







Limit cycle radius 

0.83409 

1.17243 

1.41678 

1.60137 

1.74944 

1.87276 

Thtest = 1.01 *Thcrit 







Limit cycle radius 

1.18367 

1.67068 

2.02273 

2.28975 

2.50421 

2.68354 

Thtest = 1.03 *Thcrit 







Limit cycle radius 

2.07986 

2.98424 

3.64447 

4.15197 

4.56447 

4.91321 




Fig. 8. Phase portraits for Th = 1.005Tc and various dissipative coefficient £ p . 


Qp 



c lp 




c tp 


Fig. 10. Phase portraits for Th = 1.03Tc and various dissipative coefficient £ p . 


output power. However, larger amplitude for the expansion membrane associated to higher frequencies means higher flow 
rate through the heat exchangers and regenerator. Hence, the pressure losses in addition to viscous dissipation at the 
expansion membrane reduce the efficiency. 
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Table 3 

Evolution of the performances with respect to Th and £ p . 


^P 

0.04 

0.06 

0.08 

0.1 

0.12 

0.14 

Tcriti/Tref 

1.10 

1.68 

2.09 

2.41 

2.68 

2.91 

Thtest = 1.005 *Thcrit 

Frequency (Hz) 

225 

259 

287 

309 

327 

342 

Efficiency (%) 

- 

- 

29.25 

30.79 

31.68 

32.25 

Relative efficiency 

- 

- 

55.85 

52.41 

50.37 

48.98 

Output power (mW) 

-27 

-36 

2 

55 

113 

171 

Thtest = 1.01 *Thcrit 

Frequency (Hz) 

250 

303 

345 

378 

405 

427 

Efficiency (%) 

- 

- 

29.32 

30.84 

31.71 

32.28 

Relative efficiency 

- 

- 

55.71 

52.30 

50.28 

48.91 

Output power (mW) 

-61 

-80 

17 

155 

310 

470 

Thtest = 1.03 *Thcrit 

Frequency (Hz) 

314 

421 

502 

566 

618 

662 

Efficiency (%) 

- 

- 

29.55 

31.02 

31.86 

32.40 

Relative efficiency 

- 

- 

55.19 

51.9 

49.94 

48.61 

Output power (mW) 

-237 

-253 

303 

1115 

2056 

3063 



Fig. 11. Evolution of the relative efficiency-power for various upper temperature and dissipative coefficient. 


7. Conclusion 

A nonlinear model for the analysis of the operation of a microminiature Stirling engine has been developed. 
For the understanding and the design of the engine, the nonlinear behaviour of the membrane as well as the 
thermodynamic of the engine has been considered. The existence of a Hopf bifurcation conditions the operation and 
performances. Moreover, this paper presents the centre manifold-normal form approach in order to obtain analytical 
conditions and expression for the limit cycle amplitude. This approach simplifies the analysis and the design by reducing 
the order of the dynamical system, while retaining the essential features of the dynamic behaviour near the Hop 
bifurcation point. Tough it is based on a highly idealized thermodynamic cycle the semi-analytical model is representative 
of the global behaviour of the engine and thus provides an efficient tool for preliminary design in comparison with 
extensive parametric design studies. It outlines the effect of control parameters on swept volume through the amplitude of 
the membrane flexions. Therefore, ideal performances can be evaluated and optimized. 


Appendix A 

AA. Semi-analytical method for the compression membrane 

The strategy to obtain the analytical models of the membranes is introduced with the simple example of the 
compression membrane. The geometry of the compression membrane is described in Fig. 12. 
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Fig. 12. Definition of the compression membrane parameters. 


The equilibrium equations are obtained using the Hamilton’s principle: 

(5 [ t2 (E c -E d -W ext )dt = 0 (A.l) 

where E c is the kinetic energy, E d the strain energy and W ext the work of the external forces related to the pressure in the 
considered case. H and t 2 are times such as no deflexion of the membrane occurs. 

Denoting w*(r*) for the out of plane midplane displacement of the membrane located at the radius r*, E d and W ext can be 
expressed as follows: 


E c 



w *2 r * £j r * 


(A.2) 



>ra 


7 zD 


w 


*2 

r*r 


*2 


+ 


1 

-*2 ' 


+ 


h 


2 W** 4 


r* df 


(A.3) 


where D = E/i 3 /12(1 - v 2 ) is the bending stiffness of the plate, and E, v and p are the Young’s modulus, the Poisson ratio 
and the density of the plate material, whereas h the thickness of the membrane, w** is the partial derivative to the first 
order of a function w* with respect to the variable r* and w** r * is the partial derivative to the second order of a function w* 
with respect to the variable r*. 

The total strain energy, E d of the circular plate is given as the sum of the strain energy due to bending and the membrane 
strain energy induced by large deflections. In the case of axisymmetric vibrations, the bending strain energy is given by 


■ra 


E d = E, + E nl = nD 


Wpr r * 


+ 


1 *9 


r* dr * + 


3nD 

IF 


■ra 


w*tr* 


dr" 


In the considered application, the work of the external pressure forces is 


(A.4) 


rra 

W ext = 2np / w*r*dr* (A.5) 

Jo 

where p stands for an uniform pressure applied on the surfaces of the membrane. In the case of the compression membrane 
and with respect to the notation of Fig. 1, p = ps-pc. 

If the space and time functions are supposed to be separable, the transverse displacement can be expanded in the form 
of finite series of n basis functions W*(r). Hence, with the usual summation convention for repeated indices is used over the 
range [1, n], we set w*(r,t) = g,(t) W-Xr). 

The discretization of the total strain, kinetic energy and force work expressions is made by substituting the expression 
for w*(r,t) in Eqs. (A.2)-(A.5). This leads to the following expressions: 

£( = hQj k fj 

Enl = zQiQjQkQlbyki 
E c = faqjmfj 

Wext = PQiff (A.6) 

where k~, and are the stiffness tensor, the nonlinearity stiffness tensor, the mass tensor and the force tensor, 
respectively, given by 



d 2 Wfd 2 Wt 1 dWfdWf 

L L j L L 

d r *2 d r *2 r *2 d r * dr 


r* dr* 



3D 

¥ 


■ra 


d W\ dWf 
dr* dr* 2 


+ 


1 dW\dWf 
r* 2 dr* dr 




WfWJr* dr* 


fi = 2 n 


•ra 


W'-r* dr* 


r* dr* 
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In which W*(r) = A/ 0 (/fr*) + BY 0 (/fr*) + C/ 0 (/?r*) + DJ<o(/?r*) where J 0 (x), / 0 (x) and Y 0 (x), I< 0 {x) stands for the Bessel 
functions and modified Bessel function, respectively, of the first and the second kind and p 4 = co 2 rh/D. 

Replacing E c , E d and W ext by their discretized expressions in the energy condition Eq. (A.l), calculating the derivatives 
with respect to the q{ s, and taking into account the properties of symmetry of the tensors leads to the following: 

Qi^in Qi^in ^QiQjQk^ijkn = fn (A*7) 

Eq. (A.7) represent a set of n nonlinear algebraic equations relating the n time functions q t . 

For a single vector basis we obtain 

+<?l fe ll+ 2( 7i b llll =P/l (A.8) 

where and q 1 are the acceleration and displacement response of the DOF, respectively. 

Consequently, Eq. (A.8) can be seen as a conservative Duffing’s equation. 


A.2. Semi-analytical method for the expansion membrane 


The expansion membrane is a different case because of its architecture. It is composed of two circular plates linked by a 
central rigid circular stem of radius b. The plates have two different radius a i and a 2 , respectively (Fig. 13). 

The kinetic and the total strain energy of the expansion membrane take into account the geometry. In order to simplify 
the presentation, the material and thicknesses of the two circular plates are supposed to be the same. Nevertheless, the 
methodology can be applied to a more general case. 

The transverse displacement of the lower plate which radius is a a is denoted wl*(r, t) whereas the transverse 
displacement of the upper plate is w2*(r, t). 

In the case of axisymmetric vibrations, the kinetic and the bending strain energy are given by 


fa\ 

2 1 2 3 *4’ 

ra2 

.2 1 J 2 3 ., 4 " 

E d = nD 

Jb 

W1 r** /"* 1 9 VV 1 1 -j W1 

r* 2 h 

r* dr* + nD 

Jb 

w2** r * H—yw2l +^-w2** 
r* 2 h 


r* dr* 


(A.9) 


f a 1 2 f a2 2 1 2 2 

E c = Tcph j w 1* r*dr* + 7ip/i J w2* r* dr* +-m c (wl* fa + w2* b 


The mechanical works done by the pressure forces are 

>al 


Wlext = —2npHC 


/ wl *r* dr* + - b 2 wl * h 
Jb 2 1 


= -pHCfV 


(A.10) 


(A.l 1) 


W2 ex t = 2npHC 


f® 2 1 9 

J b w2*r*dr* +-fi 2 w2* b 


= pICf2* d 


W3 ex t = 


ra\ 

2rcPatm / wVr* dr* = P atm f3* d 
Jb 


ra2 

W4 ext = -2 t ip atm / w2*r* dr* = - Patm f4* d 

Jb 


(A.12) 


(A.13) 


(A.14) 


Using the previous approach, the transverse displacements can be expanded in the form of finite series of n basis 
functions W*(r). Hence, with the usual summation convention for repeated indices is used over the range [1, n], we set 
w*(r, t) = qi(t)Wi(r). 

In the case of the expansion membrane, the local equilibrium equation for each of the circular plates is the same as for 
the compression membrane, though the boundary conditions are different. 

For the external edges we set 


Wl*,(a 1 ) = 0, W2**(a 2 ) = 0, W2**(a 2 ) = 0, W2*(a 2 ) = 0 


(A.15) 



Fig. 13. Definition of the expansion membrane parameters. 
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Table A1 


Frequency of mode 1 (Hz) 

100 

40 

20 

10 

5 

Finite element analysis 

155.401 

154.724 

149.128 

117.543 

61.196 

Analytical approach 

155.397 

154.666 

148.886 

117.051 

60.986 

Discrepancy (%) 

0.002 

0.037 

0.162 

0.418 

0.343 


Two other types of condition are related to the flexural moment and the zero slope at the central stem Wz** with z = 1, 2. 


2nbD 


wi*, r * r » +Twi* 


* r « -T-wi* 
r* 2 j | b 


+ 2nbD 


W2*« r , r « + ^W2*» r * -T-W2*, 


1 

j-* 


J lb 


= mcWl* b 


(A.16) 


Wl**(b) = 0 


(A.16a) 


W2Ub) = 0 


(A.16b) 


W\*(b) = W2*(b) 


(A.16c) 


Fig. 14 shows a scheme of the first linear flexion mode for the expansion membrane. The central rigid stem is represented 
by two punctual masses. 

In Table A1 the results form a finite element analysis and the analytical approach for al /b ratios varying from 5 to 100 
demonstrate a perfect correlation which allow the validation of the analytical approach for the evaluation of the spectral 
basis. 

Using the previously described approach for the discretized expressions of the energy terms in the Hamilton’s principle, 
the set of nonlinear algebraic equations for the expansion membrane can be obtained in the same way. 
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